#include "statistics.h"//required for statistic
#include "ClusterAnalysis.h"//required for clustering
#include "nan.h"//required for interpolation

extern "C"
{
	#include "nn.h" //required for interpolation (codes from https://code.google.com/p/nn-c/)
} 

#include <float.h>
#include <ctime>

extern const int DIME_NUM;
const double Eps = 1.0;
/** Estimate the background of a microcopy image
 *
 * @param data		Pointer to input image data
 * @param sizeX		Input image size in X direction
 * @param sizeY		Input image size in y direction
 * @param outData	Pointer to array of output image data
 * @param bitDepth	Bit depth of the input image(default 8)
 * @param binSize	Size of a tile(default 40)
 */
void backGroundEsimation(const unsigned char* data,int sizeX, int sizeY, unsigned short** outData, int bitDepth = 8,int binSize = 40)
{
	int			sizeGridX = floor((float)((sizeX-binSize/2)/(binSize/2)));
	int			sizeGridY = floor((float)((sizeY-binSize/2)/(binSize/2)));
	int			sizeGrid = sizeGridX*sizeGridY;	//number of tiles
	int*		gridX = new int [sizeGrid];		//the X coordinates of the tiling centers
	int*		gridY = new int [sizeGrid];		//the Y coordinates of the tiling centers
	float*		valueZ = new float [sizeGrid];	//the mean intensity of a tile
	double*		currentTile = new double[binSize*binSize];
	double		featureMax[DIME_NUM]={-DBL_MAX,-DBL_MAX,-DBL_MAX,-DBL_MAX,-DBL_MAX};
	double		featureMin[DIME_NUM]={DBL_MAX,DBL_MAX,DBL_MAX,DBL_MAX,DBL_MAX};
	clock_t tp1,tp2,tp3;
	vector<DataPoint> dataSets;//the array of the features
	
	/*****tiling*****/
	tp1 = clock();
	for(int i = 0;i<sizeGridY;i++)
	{
		for(int j = 0;j<sizeGridX;j++)
		{
			gridX[i*sizeGridX+j] = binSize/2*j+ binSize/2;//half binsize overlapping
			gridY[i*sizeGridX+j] = binSize/2*i+ binSize/2;
			
			for(int m = 0;m<binSize;m++)
			{
				for(int n = 0;n<binSize;n++)
				{
					currentTile[m*binSize+n] = (double)data[ (gridY[i*sizeGridX+j]-binSize/2+m)*sizeX + gridX[i*sizeGridX+j]-binSize/2+n];
				}
			}
			Statistics tileStats(currentTile,binSize*binSize);
            tileStats.calculateStatistics();
			double max = *std::max_element(currentTile,currentTile+binSize*binSize);
			double min = *std::min_element(currentTile,currentTile+binSize*binSize);
			double mean = tileStats.getMean();
			double var = tileStats.getVariance();
			double std = sqrt(var);
			double skewness = tileStats.getSkewness();
			double kurtosis = tileStats.getKurtosis();

			valueZ[i*sizeGridX+j] = mean;
			DataPoint tempDP;
			double tempPointData[DIME_NUM]={0};
			tempPointData[0] = std;
			tempPointData[1] = skewness;
			tempPointData[2] = max/min;
			tempPointData[3] = kurtosis;
			tempPointData[4] = std/mean;
			for(int k = 0;k<DIME_NUM;k++)
			{
				if(tempPointData[k] > featureMax[k]) featureMax[k] = tempPointData[k];
				if(tempPointData[k] < featureMin[k]) featureMin[k] = tempPointData[k];
			}
			tempDP.SetData(tempPointData);
			tempDP.SetDpId(i*sizeGridX+j);
			tempDP.SetVisited(false);
			tempDP.SetClusterId(-1);
			dataSets.push_back(tempDP);
		}
	}
	fprintf(stderr,"Tiling completed.\n");
	fprintf(stderr,"Use %d data points to do clustering.\n",dataSets.size());
	
	/*****clustering*****/
	
	//compute the epsilon 
	double eps=1;
	for(int k =0;k<DIME_NUM;k++)	eps = eps*(featureMax[k] - featureMin[k]);
	eps = eps * ((double)DIME_NUM+1) * 3.3234 /( (double)sizeGrid * sqrt(pow(3.1415926,DIME_NUM)) ); //3.3234 = gamma( DIME_NUM/2 +1) = gamma(3.5)
	eps = pow(eps,1/(double)DIME_NUM);
	fprintf(stderr,"Eps: %f \n",eps);

	//DBSCAN clustering
	ClusterAnalysis tileClustering;
	tileClustering.Init(dataSets,eps,DIME_NUM+1);//minPts=DIME_NUM+1.
	tileClustering.DBSCANClustering();
	tileClustering.getClusteredData(dataSets);
	vector<vector<int>> clusterPtsIds;
	tileClustering.getClusteredPointsNumber(clusterPtsIds);
	tileClustering.~ClusterAnalysis();

	//find the cluster ID with more than 200 points and minimal standard deviation
	int chooseClusterId = 0;
	if (clusterPtsIds.size()==1)
	{
		fprintf(stderr,"Only 1 class. It has %d points.\n",clusterPtsIds[0].size());
		chooseClusterId = 0;
	}
	double minstd = DBL_MAX;
	for(int i = 0;i<clusterPtsIds.size();i++)
	{
		if(clusterPtsIds[i].size()>200)
		{
			double stdmean = 0;
			for(int j=0;j<clusterPtsIds[i].size();j++)
			{
				double*  data = dataSets[clusterPtsIds[i][j]].GetData();
				stdmean = stdmean + data[0]/(double)clusterPtsIds[i].size();
			}
			if( stdmean < minstd) 
			{
					chooseClusterId = i;
					minstd = stdmean;
			}
		}
	}
	tp2 = clock();
	fprintf(stderr,"Use %d points to interpolate.\n",clusterPtsIds[chooseClusterId].size());

	/*****interpolation*****/
	int nin = 0;
	point* pin = NULL;
	int nout = 0;
	point* pinterp = NULL;
	int pcount = 0;
	point* pout = NULL;
    
	//The inputs are the grid points
	nin = clusterPtsIds[chooseClusterId].size();
	pin = (point*)malloc(nin * sizeof(point));
	for(int i = 0;i < nin;i++)
	{
		pin[i].x = gridX[clusterPtsIds[chooseClusterId][i]];
		pin[i].y = gridY[clusterPtsIds[chooseClusterId][i]];
		pin[i].z = valueZ[clusterPtsIds[chooseClusterId][i]];
	}

	//The output is the part inside the grid
	nout = (gridY[sizeGrid-1]-gridY[0]+1)*(gridX[sizeGrid-1]-gridX[0]+1);
	pinterp = (point*)malloc( nout * sizeof(point));
	for(int i = gridY[0];i < gridY[sizeGrid-1]+1;i++)
	{
		for(int j = gridX[0];j < gridX[sizeGrid-1]+1;j++)
		{	
			pinterp[pcount].x = j;
			pinterp[pcount].y = i;
			pinterp[pcount].z = NaN;
			pcount++;
		}
	}

	//interpolation
	nnpi_interpolate_points(nin, pin, -DBL_MAX, nout, pinterp);

	//copy the output for extrapolation
	pout = (point*)malloc(sizeX*sizeY * sizeof(point));
	for(int i = 0;i < sizeY;i++)
	{
		for(int j = 0;j < sizeX;j++)
		{
			pout[i*sizeX+j].x = j;
			pout[i*sizeX+j].y = i;
			pout[i*sizeX+j].z = NaN;
		}
	}
	nout = 0;
	for(int i = gridY[0];i < gridY[sizeGrid-1]+1;i++)
	{
		for(int j = gridX[0];j < gridX[sizeGrid-1]+1;j++)
		{	
			pout[i*sizeX+j].z = pinterp[nout].z;
			nout++;
		}
	}

	/*****extrapolate*****/
	for(int i = gridY[0];i < gridY[sizeGrid-1]+1;i++) //every line
	{
		int head = -1;	//index of the first none NaN point
		int end = -1;	//index of the last none NaN point
		for(int j = 0;j < sizeX-1;j++)
		{
			if( _isnan(pout[i*sizeX+j].z) && !(_isnan(pout[i*sizeX+j+1].z))) 
			{
				if (end-head<sizeX/2) head = j+1;
			}
			if(!(_isnan(pout[i*sizeX+j].z)) && _isnan(pout[i*sizeX+j+1].z) )
			{
				if (head > end)	end = j;
			}
		}
		if (end-head < sizeX/2) continue;	//do not extrapolate when the useful part is too short

		double k1 = pout[i*sizeX+head+1].z - pout[i*sizeX+head].z;
		double k2 = pout[i*sizeX+end].z - pout[i*sizeX+end-1].z;
		for(int j = 0;j < head;j++)
		{
			pout[i*sizeX+j].z = pout[i*sizeX+head].z - k1 * (head - j);
		}
		for(int j = end +1 ;j < sizeX;j++)
		{
			pout[i*sizeX+j].z =pout[i*sizeX+end].z + k2 * (j - end);
		}
	}
	for(int j = 0;j < sizeX;j++)		//every colunm
	{
		int head = -1;
		int end = -1;
		for(int i = 0; i<sizeY-1;i++)
		{
			if(_isnan(pout[i*sizeX+j].z) && !(_isnan(pout[(i+1)*sizeX+j].z))) 
			{
				if(end-head <sizeY/2) head = i+1;
			}
			if(!(_isnan(pout[i*sizeX+j].z)) && _isnan(pout[(i+1)*sizeX+j].z)) 
			{
				if(end < head)end = i;
			}
		}
		if (end-head < sizeY/2) continue;

		double k1 = pout[(head+1)*sizeX+j].z - pout[head*sizeX+j].z;
		double k2 = pout[end*sizeX+j].z - pout[(end-1)*sizeX+j].z;

		for(int i = 0;i < head;i++)
		{
			pout[i*sizeX+j].z = pout[head*sizeX+j].z - k1 *(head-i);
		}
		for(int i = end+1;i < sizeY;i++)
		{
			pout[i*sizeX+j].z = pout[end*sizeX+j].z + k2*(i - end);
		}
	}
	tp3 = clock();

	/*
	FILE *logfile;
	logfile = fopen("C:\\Users\\Lm\\Documents\\MATLAB\\Background estimation\\logfile.txt","a");
	fprintf(logfile,"Clustering running time is %f \n",(double)(tp2 - tp1) * 1000.0 / CLOCKS_PER_SEC);
	fprintf(logfile,"Interpolation running time is %f \n",(double)(tp3 - tp2) * 1000.0 / CLOCKS_PER_SEC);
	fclose(logfile);//*/
	cout<<"Clustering running time is "<<(double)(tp2 - tp1) * 1000.0 / CLOCKS_PER_SEC << " ms!" << endl;
	cout<<"Interpolation running time is "<<(double)(tp3 - tp2) * 1000.0 / CLOCKS_PER_SEC << " ms!" << endl;//*/
	
	//check the results
	int nacount = 0;
	unsigned char imageMin = *std::min_element(data,data+sizeX*sizeY);
	unsigned char imageMax = *std::max_element(data,data+sizeX*sizeY);
	for(int i = 0;i < sizeY;i++)
	{
		for(int j = 0;j < sizeX;j++)
		{
			if (_isnan(pout[i*sizeX+j].z)) 
			{
				nacount++;
				pout[i*sizeX+j].z = 0; //set it to zero if NaN
			}

			//The background should not be larger than maximum.
			if(pout[i*sizeX+j].z > imageMax)
				pout[i*sizeX+j].z = (double)imageMax;

			//Should not be smaller than the minimum.
			if(pout[i*sizeX+j].z < imageMin) 
				pout[i*sizeX+j].z = (double)imageMin;
		}
	}
	fprintf(stderr,"%d points NaN.\n",nacount);

	*outData = (unsigned short*)malloc(sizeX*sizeY*sizeof(unsigned short));
	for(int i = 0;i<sizeX*sizeY;i++)
	{
		unsigned short * data;
		data = &(*outData)[i];
		if (bitDepth == 8)
			*data = (unsigned short)(256 * pout[i].z);//output 16bit data
		else if (bitDepth == 16) 
			*data = (unsigned short)pout[i].z;
	}
	delete	gridX;
	delete	gridY;
	delete	valueZ;
	delete	currentTile;
	free(pin);
	free(pinterp);
	free(pout);
}